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Abstract. Wc analyse an eleven-orbit inspiral of a non-spinning black- hole 
binary with mass ratio q = M1/M2 = 4. The numerically obtained gravitational 
waveforms are compared with post-Newtonian (PN) predictions including several 
sub-dominant multipolcs up to multipolar indices (I = 5, m = 5). We find that 
(i) numerical and post-Newtonian predictions of the phase of the (2, 2) mode 
accumulate a phase difference of about 0.35 rad at the PN cut off frequency 
Mu> = 0.1 for the Taylor Tl approximant when numerical and PN waveforms 
are matched over a window in the early inspiral phase; (ii) in contrast to 
previous studies of equal-mass and specific spinning binaries, we find the Taylor 
T4 approximant to agree less well with numerical results, provided the latter 
are extrapolated to infinite extraction radius; (iii) extrapolation of gravitational 
waveforms to infinite extraction radius is particularly important for subdominant 
multipolcs with I 7^ m; (iv) 3PN terms in post-Newtonian multipole expansions 
significantly improve the agreement with numerical predictions for sub-dominant 
multipolcs. 

The research area of Gravitational Wave (GW) Physics has acquired enormous 
momentum in the course of the last decade. The ground-based detectors LIGO, 
VIRGO and GEO600 have operated at design sensitivity and the former two are 
currently being upgraded to advanced status [H [2]. Simultaneously the planned 
ES A/NASA space mission LISA will soon enter a crucial stage with the launch 
of the precursor mission LISA Pathfinder [3J. Progress in theoretical GW source 
modelling has mirrored that on the experimental side. Most importantly, the general 
relativistic two-body problem for comparable-mass systems has been solved using 
numerical relativity (NR) techniques [H [5j [6] , leading to a wealth of insight into the 
dynamics of black- hole binaries [3 [8] . 

In order to maximise the scientific output of the expected gravitational wave 
observations, it is crucial to have available catalogues of highly accurate gravitational 
waveform templates which are then used via the matched filtering technique [3] to dig 
out physical signals from the observed data stream. The construction of such waveform 
templates currently follows either of the two following strategies. Phenomenological 
template banks are based on hybrid waveforms combining post-Newtonian predictions 
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for most of the inspiral with numerical relativity results for late inspiral, merger 
and ringdown. Phase and amplitude are then assumed to be well-approximated by 
relatively simple functions of the binary parameters and the GW frequency u. Free 
parameters in these functions are determined by comparison with a finite number of 
hybrid waveforms. This approach, initially presented in [lOj [TTJ [12] has been applied 
to a subset of spinning configurations in [13j HU [HI [16]. The second approach, 
the effective-one-body method (EOB) [17j . targets at semi-analytic predictions for 
the entire waveform via matching post-Newtonian results to an effective-one-body 
metric and thus models the binary through merger and ring-down. Again, free 
parameters in the matching arc determined by comparison with a finite set of numerical 
relativity simulations. Applications of the EOB method can be found, for example, in 
[T51 IT91 [201 |2"TI 122] for non-spinning and in [33] for non-precessing spinning binaries. 

Clearly, both approaches are heavily dependent on high-precision numerical 
relativity results. Length requirements for the numerical simulations in the context of 
gravitational wave detection have been studied for non-spinning binaries and special 
spin configurations in [24] and are in the range of about ten orbits. A recent 
study on accuracy requirements for detection efficiency and parameter estimation 
by MacDonald et al. [25] employs criteria developed in [26] [22] and reports a larger 
number of about 30 orbits required for hybridisation with current post-Newtonian 
waveforms; for more details see these articles and references therein. Finally, we note 
that numerical relativity waveforms are directly used in GW data analysis inside the 
Ninja project [27 ] [28 ] , 

Purpose of our paper is to study in detail the accuracy of an eleven-orbit inspiral 
of a non-spinning black-hole binary with mass ratio q = 4 obtained with the Lean 
code [29 . While higher-order multipoles have been investigated in the context of the 
EOB model [30j [31] and PN-NR comparisons in [3"2] [33 ] [33 ] [35], we are not aware 
of their inclusion in the construction of phenomenological models of hybrid PN-NR 
waveforms. In this work, we will compare our numerical results for the quadrupole 
as well as several subdominant modes with the Taylor Tl and T4 approximants and 
different PN multipolar expansions. 

We start this paper with a brief summary of the numerical framework in Sec. [TJ 
The simulations and the numerical error analysis is presented in Sec. [2] In Sec. [3] we 
discuss the different post-Newtonian approximants used in this work and the method 
to construct hybrid waveforms. These are analysed in Sec. [Hand we conclude in Sec. [5] 

1. Computational framework 

The initial data for our black-hole binary configurations are constructed according 
to the puncture method [36]. Specifically, we use the analytic Bowen-York extrinsic 
curvature [37. with linear momentum P and solve the Hamiltonian constraint for the 
conformal factor using the spectral solver of Ansorg et al. 38J . In order to reduce the 
eccentricity of the binary system, we determine the non-vanishing radial component 
-Prad of the initial momentum via the iterative procedure described in Ref. [39j . 




The initial parameters thus obtained are given in units of the total black-hole 
mass M = Mi + in Table [T] The bare mass parameters have been chosen such 
that the irreducible masses as calculated with Thornburg's AHFinderDirect [40). [41] 
correspond exactly to a binary with mass ratio q = 4. For completeness we also give 
the final spin of the merged hole and the recoil velocity. The latter is in excellent 
agreement with the prediction t>kick = 156.9 km/s by [42], confirming that the early 
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Table 1. Initial bare mass parameters mi, rri2, location asi, £2 and tangential 
and radial linear momentum of black holes 1 and 2 of the binary, v-^^ and ja n 



are the kick velocity and spin of the post-merger hole. 


mi/M m 2 /M xi/M x 2 /M P taa /M P ia d/M 


Wkick 


jfin 


0.7923 0.1913 2.1865 -8.746 ±0.05805 T3.894 x 10~ 4 


156.6 km/s 


0.473 



inspiral does not significantly contribute to the recoil; cf. the discussion in Sec. Ill C 
of [43]. 

These initial data are evolved in time with the so-called moving puncture 
framework [5] [3] using the Lean code [2H] with upgrades to sixth-order spatial 
discretization as summarised in Sec. Ill of Ref. [44]. The Lean code is based on the 
Cactus computational toolkit [45] and employs mesh refinement provided by Carpet 
[46] [47] . The exact implementation of the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) [48] [49] formulation of the Einstein equations is given by Eqs. (Al), (A4), 
(A6-A8) of [25] ■ We further impose a floor value Xfloor = 10~ 4 on the conformal 
factor; cf. [5]. The gauge variables are initialised as vanishing shift f3 l = and a 
precollapsed lapse a = y/x, where x = (det7y )~ 1 / 3 and 7^ is the three metric. Lapse 
and shift are evolved according to Eqs. (17) and (26) of Ref. [50] with a damping 
parameter r\ = 1.75/M. 

The computational domain consists of a set of five nested boxes centred on the 
coordinate origin and five additional refinement levels with two components each, 
centred on the black holes. In terms of the notation in Sec. II E of [29], the exact grid 
setup is given in units of M by 

{(307.2, 153.6, 102.4, 32, 16) x (3.2, 1.6, 0.8, 0.4, 0.2), hi} 

with resolutions h Q = Af/160, hi = M/180, h 2 = M/200, h 3 = A//220 and 
hi — Af/240 for studying the convergence properties. The lowest resolution simulation 
using ho becomes unstable due to a gauge instability at the outermost refinement 
boundary at about t — 1700 M; see [5T] [521 [S3] 03] for more in-depth discussion and 
cures of this instability. While we do not use this simulation directly for comparison 
with post-Newtonian predictions, it extends sufficiently far into the inspiral to provide 
consistency checks of our convergence results. 

Gravitational waves are extracted in the form of the Newman-Penrose scalar ^4 
at radii r cx = 64 M, 72 M, 80 M, 88 M, 96 M, 104 M and 112 M using the method 
described in Appendix C in [29) . The resulting ^4 is decomposed into multipoles ipi m 
by projection onto spherical harmonics of spin- weight —2 according to Eq. (2.1) of 
Ref. [33]; see also [55] for sign conventions. 

Our analysis of the gravitational waves will use the gravitational wave strain h 
which is related to the Newman-Penrose scalar by 

y i = h + -ih x . (1) 

The multipoles Hi m of the strain are related to those of ^4 by Eq. (II. 5) of Brown et 
al. [56], that is, two integrations in time of ^i m - These integrations in time represent 
a non-trivial operation and often result in low frequency drifts [331 157) . In order to 
circumvent these problems, we use a modified version of the integration in Fourier 
space as originally introduced in Ref. [58], see also [35]: the Fourier transform Hi m {u) 
is divided by — w 2 , but set to zero inside a specified window. Our modification consists 
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Figure 1. Convergence of the phase tf>22 for resolutions hi, h,2, /13 and /14 and 
aligning the waveforms at the maximum of the (2, 2) mode (left panel) and for 
resolutions ho, hi, h^, ho, and /14 and aligning the waveforms at the start of the 
numerical simulation (right panel) . The scaling factors correspond to eighth-order 
convergence. The lower solid curve in the left panel represents the amplitude A22 
and serves as orientation. 



in smoothing the filter function from Heaviside shape to 

[ x < 

9( x ) = \ W I x4 ( 1 ~ x ) 4dx < x < 1 ( 2 ) 

[l x > 1 

where x = (u> — ljq)/ Aoj and AT is a normalisation constant. For the present study, we 
empirically find AIljq = 0.001 and M Au> = 0.001 to provide satisfactory results. 

2. Numerical results 

Before analysing the physical properties of the binary evolution, we estimate the 
uncertainties due to discretization and wave extraction at finite radii. The convergence 
of the phase <f> obtained from the I = 2, m = 2 multipole is shown in Fig. Q] 
for all resolutions hi, i = 0,...,4. For the results in the left panel, we have 
aligned the waveforms in time at tA 22 , the time of the peak of the amplitude of 
the (2,2) multipole. Those displayed in the right panel are obtained for aligning 
the waveforms at to, the start of the simulatiory]. The results are consistent with 
eighth-order convergence for all resolutions employed. Similarly, we observe eighth- 
order convergence for the amplitude. Bearing in mind that most of the discretization 
in our code is sixth-order, it appears that the leading-order discretization errors are 
subdominant. Whether this occurs due to systematic cancellation is hard to identify 
because of the enormous complexity of the Einstein equations. In order to test whether 
the observed convergence is coincidental or systematic, we have performed the larger 
set of runs employing five different resolutions and using different choices of aligning 
the waveforms in time. The observation of the same clean convergence order for 
all cases demonstrates the robustness of our observations. For the determination of 
uncertainties due to discretization, however, we will calculate Richardson extrapolated 
results using a more conservative sixth-order extrapolation. 

The deviations of phase and amplitude obtained at finite resolutions from the 
Richardson extrapolated values are shown in Fig. [2] for both types of alignment of the 



X The lowest resolution simulation has only been used in the latter case because it does not extend 
to the maximum of the amplitude of the (2,2) mode. 



Eleven-orbit inspiral of a mass ratio 4-tl black-hole binary 5 




(t-RJ/M (t-RJ/M 

Figure 2. Deviation of the phase <f> (left) and relative deviation of the amplitude 
(right panel) obtained for finite resolution from the Richardson extrapolated 
values and for aligning the waveforms at *a 2 2 > ^ ne maximum of the amplitude 
of the (2,2) multipole (upper panels) and at the start of the simulation to (lower 
panels). The vertical lines line marks the time where the frequency of the (2,2) 
multipole reaches Mu>22 =0.1. 

waveforms, at tA 22 (upper panels) and at to (lower panels). For the former alignment, 
the figure implies a phase error of A(j) < 0.6 rad for the high-resolution simulation. 
We note, however, that the phase difference between the high resolution and the 
extrapolated result is nearly constant until the late ringdown stage. Because a constant 
phase offset does not enter the comparison with PN results, the relevant phase error 
for this purpose is given by the variation of the difference over the simulation which 
is significantly smaller, about A(f> ss 0.11 rad. The relative amplitude error shown 
in Fig. [2] is less than 1 % and, as we shall see below, dominated by the uncertainty 
arising from the use of finite extraction radii. 

In case of aligning the waveforms at to, Fig. [2] implies an accumulated phase error 
of about 0.4 rad for the phase of the high resolution run at i w defined as the time 
when the frequency of the multipole reaches Mu>22 = 0.1, the endpoint of the post- 
Newtonian integration in Sec. [3] (vertical lines in the figure). The amplitude error has 
grown to 2 % at the same time. We note here that alignment of the gravitational waves 
at tA 22 results in smaller uncertainty estimates. For the analysis of Sec. HI however, 
estimates obtained for both type of alignments are adequate. In that analysis we will 
exclusively use the high-resolution data set. 

In order to estimate the error due to finite extraction radius, we assume a 
polynomial dependence of the phase error on l/i? C x [59] . While the linear term is 
expected to dominate the error, in practise the quadratic term may also be significant 
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Figure 3. Deviation of the phase 4> (left) and relative deviation of the amplitude 
(right panel) obtained at finite extraction radius from the values extrapolated 
according to Eq. J3}. 



[50] . We therefore fit the numerical data assuming either of 

f{r) = k + f ~, (3) 
r 

/W = /o + 7 + ^ ( 4 ) 

and denote the extrapolated results thus obtained "xpoll" and "xpol2" , respectively. 
Eventually, we use the extrapolated result "xpoll" in the analysis below and estimate 
the error by its deviation from the "xpol2" result. These deviations are shown for 
phase and amplitude and a subset of all extraction radii in Fig. [3] From the figure we 
infer a relative error of 5 % or less for the amplitude while the phase error is < 0.2 rad 
up to merger and increases to about 0.5 rad in the late ringdown stage. 

The phase uncertainties of higher order multipoles behave similarly to those of the 
quadrupolc, but we find the overall errors to scale approximately with the multipole 
index m. The amplitude errors of subdominant modes are significantly larger, however, 
due to the lower signal combined with numerical noise. Extrapolation of the amplitude 
to infinite extraction radius amplifies the numerical noise and we therefore present 
higher order multipoles using extrapolated phase, but amplitude from the largest 
extraction radius instead. The numerical uncertainties of these amplitudes are about 
12 % for the (3,3) mode, 20 % for the (4,4) and (5,5) mode and 25 % for the (3,2) 
and (4, 3) mode. While these errors are too large for a high precision study, they 
will allow us to calibrate the significance of 3PN order terms in the post-Newtonian 
multipole expansion in Sec. HI 

The errors stated so far are internal checks of accuracy. For the purpose of 
an external verification of our results, we compare the amplitude and phase of the 
1 = 2, m = 2 mode to results obtained with the independent BAM code [BTJ [55] . In 
order to avoid additional uncertainties arising from the integration of the Newman- 
Penrose scalar $4 to strain h, here we compare the modes of ^4 which we directly 
obtain from the numerical simulations in both codes. The BAM simulations also 
use lower resolutions (hence the larger uncertainties) and different extraction radii. 
We therefore present this comparison using results extrapolated both in resolution, 
assuming second order convergence for the BAM runs, and extraction radius, using 
Eqs. ©, ®. 

The uncertainties of the BAM simulation are studied in detail in [62] and are 
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Figure 4. Relative amplitude differences (left panel) and phase differences (right 
panel) of the 1 = 2, m = 2 modes of ^4 between the extrapolated (in resolution 
and extraction radius) results from the BAM and the Lean codes. Results from 
both codes are aligned at the time of peak amplitude tj\ 22 . We note that the 
relatively large percentage discrepancy between the results of the two codes at 
late times is a consequence of the very low amplitude of the wave signal in the 
late ringdown. 



summarised as follows. The accumulated phase error for the N — 96 high resolution 
run at as obtained from Richardson extrapolation and for aligning the waveforms 
at to is 0.46 rad and the error due to the use of finite extraction radius at i? ox = 90 M 
is 0.15 rad. The corresponding uncertainty in the amplitude is about 10 % for both, 
discretization and finite extraction radius during the inspiral and merger. 

Relative deviations between the amplitudes and phase differences obtained from 
the two codes are shown in Fig. |4j Note that the waveforms have been aligned at 
£^22 defined as t = on the horizontal axes of this figure. This alignment becomes 
necessary because of the different initial separations (and, hence, duration) of the 
simulations, D — 10 M for BAM and D = 11 M for Lean. Phase and amplitude 
differences shown in the figure are well within the uncertainty estimates of the two 
simulations. 

Finally, we evaluate the residual eccentricity of the binary simulated with Lean 
as a function of time according to [63] 

,.x 0NR.(Q - <fetft) f .s 

e<t>{t) = , (5) 

where </>nr is the phase of the numerical (2, 2) mode and (f>fn a seventh order polynomial 
fit of this phase over a time window t — R cx = 250 . . . 1600 M . We thus obtain an 
oscillatory pattern of e^i) with an amplitude gradually decreasing from 0.005 early 
on towards 0.002 one to two orbits before merger. 

3. Post-Newtonian approximants and hybrid waveforms 

Post-Newtonian waveforms are calculated using the so-called Taylor Tl and Taylor 
T4 approximants; cf. |59j for a summary including further approximants. Frequency 
and phase are obtained from the system of ordinary differential equations 

dx L , . 
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Here x is related to the orbital frequency via x = (Mf2) 2 / 3 , <i> is the orbital phase, 
and the gravitational wave flux C and the orbital energy E are given by Eqs. (231) 
and (203), respectively, of Blanchet's review [64,. The difference between the Taylor 
Tl and Taylor T4 approximant is the treatment of the right hand side of Eq. ([5]). 
The former calculates the quotient C/(dE/dx) numerically from the PN truncated 
individual expressions for numerator and denominator while the T4 approximant 
expands the quotient in a series in x and truncates this series at the appropriate 
PN order, 3.5 in our case. The two approximants thus differ only at higher PN order, 
but previous studies have found the Taylor T4 approximant to produce particularly 
good agreement with numerical results [59j [62] . 

In either case, we obtain the rescaled frequency x and the orbital phase $ from 
integration of Eqs. ^ and ([7]). The gravitational wave multipoles are given in terms 
of these quantities by Kidder's Eqs. (79)-(116) in [55] (to 3PN order for (2,2) and 
(4,4) and 2.5PN order otherwise) or Eq. (9.4) of the more recent study by Blanchet 
et al. 66J (to 3PN order). Unless specified otherwise, results displayed use the latter 
3PN terms. 

In order to construct a hybrid waveform, we need to determine the integration 
constants <&o and to of the system ([6]) , (0 . In practise, this is achieved by maximising 
the overlap between the real part of the post-Newtonian multipole 7?22,pn and its 
numerical counterpart i?22,NR over a matching window t\ < t — R cx < ti M . In this 
context, the overlap of two functions / and g is defined as 

(f,g)= ( 2 f(t)g(t)dt. (9) 



The maximisation is achieved using the downhill simplex method of Nelder & Mead 
|67[ 168] . Finally, we combine the PN and numerical waveform into a hybrid according 
to 

Him = (1 - w)<xff/ m ,PN + wH hn ^ R , (10) 

where the weighting function w = 1 for t < ti, w = for t > £2 and for values 
ti < t < ti a smooth transition is given by 

w=630 (l z9 -r s+ r 7 -r 6+ r 5 )' 



t - 1 



(i2 > 

The additional factor a has been introduced to compensate for differences in the 
amplitude between the post-Newtonian and numerical results. For each multipole it 
is chosen such that the average PN wave amplitude inside the interval [ti, £2] matches 
the NR result. 

In order to test the robustness of our results versus the details of the matching 
procedure, we perform an alternative matching by equating phase and time at a 
fiducial point in time chosen to be t^ where the gravitational wave frequency of the 
(2,2) mode takes the value Mw 2 2 = 0.1. 
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Figure 5. Real part of the 1 = 2, m = 2 and the 1 = 3, m = 3 gravitational wave 
modes. The Taylor Tl PN waveforms (dashed) are matched to the numerical 
results in the window 250 M < t — i? cx < 600 M indicated by vertical dotted 
lines. 

4. Analysis 

The results of the matching procedure using a window t\ — 250 M, t 2 = 600 M are 
illustrated for the Taylor Tl approximant in Fig. [5j which shows the real part of the 
(2, 2) and the (3, 3) modes. The matching window is indicated by the vertical dotted 
lines. The corresponding figure for Taylor T4 would look nearly identical. 

In order to quantitatively analyse the agreement between the PN and NR results, 
we display in Fig. [5] the phase differences for various multipoles obtained by using 
the Taylor Tl (left panels) and T4 (right panels) approximants and numerical results 
extracted at R cx = 64 M and 96 M as well as extrapolated to infinite extraction 
radius. For orientation we plot near the bottom of each panel the gravitational wave 
frequency Mw 2 2 of the 1 = 2, m = 2 mode (dash-dash-dotted curve). All curves end 
at Mui22 = 0.1 which is the maximum frequency chosen for the PN integration. 

We begin our discussion with the dominant quadrupole mode I = 2, m = 2. In all 
cases, the phase agreement between PN and NR results is better than Acj) = 0.1 rad 
inside the matching window, but gradually increases as the inspiral proceeds until 
the accumulated phase discrepancy reaches values between 0.25 rad and 0.75 rad. At 
both finite extraction radii the agreement grows to larger values for the Taylor Tl 
expansion than the T4 approximant. The behaviour of the (3,3), (4,4) and (5,5) 
modes is similar, although with larger phase discrepancy approximately proportional 
to the multipole index m. This is not too surprising bearing in mind that the phase 
grows faster for higher-order modes with the same proportionality. 

Those subdominant multipoles with I =^ m, however, exhibit a drastically different 
behaviour when extracted at finite radius. For both, the (3, 2) and the (4, 3) mode the 
numerical phase differs significantly from the PN predictions throughout the entire 
simulation. We note, in this context, that the orbital phase in the PN expressions 
is fixed exclusively using information from the (2, 2) mode. Furthermore, inspection 
of the post-Newtonian multipoles in Refs. [65] and [66] reveals that the (2,2) and 
the (3, 2) GW multipoles are almost in phasdjj. Indeed, this is also the case for the 

§ The minor dephasing due to the complex PN amplitudes of the (2, 2) and (3, 2) modes is negligible 
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Figure 6. Phase difference between the PN and NR results obtained by using 
a matching window with ti = 250 M and t2 = 600 M and numerical results at 
(from top to bottom) R cx = 64 M, 96 M and extrapolated to infinity. Results in 
the left column are obtained for the Taylor Tl expansion, those on the right for 
Taylor T4. The dash-dash-dotted curve near the bottom of each panel gives the 
GW frequency of the (2, 2) mode for reference. 



numerical modes provided we extrapolate to infinite extraction radius, as becomes 
apparent in the bottom panels where the large phase error of the (3, 2) and (4, 3) 
modes disappear. We also illustrate this feature in Fig. [7] where we plot the real 
part of the numerical (2,2), (3,2) and (4,2) mode^jj. The figure demonstrates that 
(i) the dephasing is largest for the (4, 2) mode, (ii) the dephasing decreases at larger 
extraction radius and (iii) the dephasing virtually disappears if we extrapolate to 
infinite extraction radius. We make exactly the same observations for the (4, 3) and the 
(3,3) multipoles; see the short-dashed and dash-dotted curves in Fig. [SI Our findings 

in this context. 

|| The low amplitude of the (4, 2) combined with numerical noise prevent us from performing accurate 
extrapolation to infinite extraction radius which is also the reason it is not included in the remainder 
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Figure 7. The (2,2) (solid), (3,2) (dashed) and (4,2) (dotted) multipole 
extracted at R e x = 64 M and 96 M (upper panels) and the former two 
extrapolated to infinite extraction radius (bottom panel). 

demonstrate the importance of having available either reliable wave extraction tools 
at future null infinity 69 , 70] 171] or accurate extrapolations from results at finite radii 

A further interesting feature of the extrapolated modes is the significant 
improvement of the phase agreement between NR and the Taylor Tl approximant; for 
R cx — > oo, Taylor Tl provides better agreement than T4. We thus cannot confirm for 
our unequal-mass binary the exceptional behaviour of the T4 approximant reported 
for the equal-mass case by [59] and a q — A unequal-mass binary in Fig. 8 of Ref . [62 . 
We note, however, that this does not constitute an incompatibility of our results with 
those of Ref. [62] , because the latter are based on waveforms at finite extraction radius 
90 M, where we also obtain better agreement of numerical results with T4. 

We next perform the following tests in order to investigate the robustness of 
our observations at infinite extraction radius, (i) Motivated by the observation of 
noise in the numerical phase due to spurious initial radiation in the earlier part of 
their waveform in Ref. [21], cf. their Fig. 17§, we choose a later matching window 
ti = 900 M, t% = 1250 M. (ii) Instead of using a window, we match phase and time 
at the fiducial point in time t u where Mu>22 =0.1 as done for example in Ref. [62 ]. 
(iii) We apply the matching procedure with t\ = 250 M, t 2 = 600 M to the BAM 
waveform. The results are shown, from top to bottom, in Fig. [8] As expected, different 
alignment procedures produce different functions of time for the phase differences. 
But for all cases, we observe better agreement of the numerical simulations with the 
Taylor Tl prediction as compared with Taylor T4. This is also the case for the BAM 
simulation which is furthermore consistent within the respective uncertainties with 
that obtained with the Lean code. We conclude that the exceptional agreement of 
Taylor T4 observed for some specific black-hole configurations is most likely explained 
by a coincidental cancellation of higher-order post-Newtonian corrections which does 
not hold for general systems. 

of this study. 

1 We note that Boyle et ctl. use different initial data and that the amplitude of the noise is at least 
one order of magnitude below the effects investigated in our work, compatible with the absence of 
such noise in the (2,2) mode on the scale of our Fig. [6] 






-°4 



600 

(t-RJ M 



Figure 8. Phase difference between the PN and NR results obtained by using 
numerical results extrapolated to infinity and employing a later matching window 
tl = 900 M , t2 = 1250 M (upper panels), matching phase and time at a fiducial 
point in time where Muiii = 0.1 (centre panels) and using a matching window 
ii = 250 M, t 2 = 600 M for (2,2) and (3,2) multipoles of the BAM simulation 
(bottom panels). Results in the left column are obtained for the Taylor Tl 
expansion, those on the right for Taylor T4. The dash-dash-dottcd curve near 
the bottom of each panel gives the GW frequency of the (2, 2) mode for reference. 



Finally, we consider differences in the amplitude predictions of post-Newtonian 
and numerical relativity results. The PN expressions for the gravitational wave 
multipoles in Refs. 65 and 66 differ in the inclusion of higher order terms in several 
subdominant modes in the latter work. In order to assess the significance of these 
higher-order terms, we have performed the matching of numerical to PN waveforms 
as described in Sec. [3] with t\ — 250 M and i2 = 600 M using either set of multipole 
expressions. As mentioned above, the matching process involves a rescaling of the PN 
waveforms to compensate for amplitude differences in the matching window, the factor 
a in Eq. (fTOj) . For all subdominant multipoles which have been extended to 3PN order 
in [66] , we observe in Table [2] an improvement in the agreement between NR and PN 
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Table 2. Deviations of the rescaling factor a introduced in Eq. II10I I from unity 
for several multipoles and using the multipole expressions of Ref. 1651 or Ref. I66| 
which includes higher-order PN terms for subdominant modes. 



(l,m) 


(2,2) 


(3,3) 


(3,2) 


(4,4) 


(4,3) 


(5,5) 


a — 


1 (Tl, [55]) 


0.047 


-0.071 


-0.065 


-0.108 


0.251 


0.712 


a — 


1 (Tl, [56]) 


0.047 


0.023 


-0.059 


-0.108 


-0.075 


-0.106 


a — 


1 (T4, [55]) 


0.048 


-0.071 


-0.065 


-0.107 


0.251 


0.712 


a — 


1 (T4, [55J) 


0.048 


0.023 


-0.058 


-0.107 


-0.075 


-0.106 



amplitudes, i. e. a value a closer to unity. Even bearing in mind the relatively large 
uncertainties in the numerically obtained amplitudes, this improvement is significant, 
at least for the (4, 3) and (5, 5) mode. 

5. Conclusions 

We have studied numerical simulations of a non-spinning black-hole binary system 
with mass-ratio q = 4 lasting about 11 orbits prior to coalescence. The numerical 
uncertainties for phase and amplitude due to discretization are A(f> ss 0.6 rad and 
AA/A r; 1 % for the quadrupole mode through inspiral, merger and ringdown when 
aligning the waveforms at the peak of the amplitude of the (2,2) multipole. The phase 
error is approximately constant, however, and we estimate the uncertainty for the 
purpose of a PN comparison closer to A(f> sa 0.11 rad. Numerical error due to wave 
extraction at finite radii results in larger uncertainties for the amplitude of about 
AA/A < 5 % and a phase error Acf> < 0.2 rad up to merger and 0.5 rad in the late 
ringdown stage. Uncertainties for subdominant multipoles are larger; approximately 
in proportion to the multipole index m for the phase and reaching 10 % to 25 % for 
higher-order modes in the amplitude. We also observe agreement within numerical 
uncertainties with an independent simulation of a q — 4 binary obtained with the 
BAM code. 

We have performed a matching to post-Newtonian predictions using the Taylor Tl 
and T4 approximants and employing 3PN accurate expressions for the GW multipoles. 
Our main results in this comparison can be summarised as follows. 

Using our numerical waveforms at finite extraction radii gives the misleading 
impression that Taylor T4 produces better agreement than the Tl approximant. 
Including subdominant multipoles with I m explicitly demonstrates the internal 
inconsistency of the numerical results at finite (too small) extraction radii: contrary 
to expectations the (2, 2) and (3, 2) as well as the (3, 3) and (4, 3) multipoles are 
significantly out of phase. This feature improves when going to larger radii and 
disappears when results are extrapolated to infinite radius. Subdominant modes thus 
provide valuable tests for the internal consistency of numerical results, irrespective of 
whether they are included in a comparison with post-Newtonian predictions or not. 

By using extrapolated numerical results, we find the Taylor Tl approximant to 
provide better agreement with numerical results: using a matching window in the 
early inspiral, the accumulated phase error at Mwii = 0.1 is A<fi « 0.35 rad for Tl 
compared with 1.0 rad for T4. For subdominant multipoles we obtain deviations of 
PN from NR results which to good approximation are proportional to the multipole 
index m. 

The inclusion of higher-order PN terms in expressions for sub-dominant 
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multipoles in 66 leads to improved amplitude agreement with numerical results in 
our matching window covering approximately the frequency range 0.05 < Mu> < 0.06. 
Even bearing in mind the relatively large numerical errors for the low amplitude modes, 
this improvement is significant at least for the (4,3) and the (5, 5) multipole. 
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